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1 Introduction 

Networks of coupled dynamical systems play an important role for all branches 
of science [H Ej El H] . In the neuroscience, for instance, there is a need for 
modeling large populations of coupled neurons in order to approach problems 
connected with the synchronization of neural cells or other types of collec- 
tive behavior [5l [6l H]. The investigation of the dynamics of coupled lasers 
[3 El El do] is important for many purposes including secure communication 
[TTl [T2] or high-power generation. The interacting biological, mechanical or 
electrical oscillators [131 [E] belong already to classical models for studying 
various aspects of collective dynamics. In neural networks, the synchronous 
activity might be pathological [15], and hence, there was recently an in- 
creasing effort to control the desynchronization of populations of coupled 
oscillators. In particular, the coordinated reset stimulation technique [H [32] 
proposes to establish a cluster-state in the network, in which the oscillator's 
phases split into several subgroups. This example illustrates the importance 
of the analysis of cluster formation in coupled systems. Our paper inves- 
tigates the connection between the properties of a single oscillator, i.e. its 
sensitivity to stimulations, and the formation of clusters in a globally cou- 
pled system of such oscillators. We show that by altering the shape of the 
sensitivity function, called the phase response function, different clusters in 
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a network can be stabilized. More presicely, we study a family of the phase 
response curves, which are unimodal and turn to zero at the spiking moment. 
This choice is motivated by several well known neuron models. It appears 
that the position of the maximum of the unimodal sensitivity function with 
respect to the spiking point plays an important role for determining whether 
the system will synchronize or approach a two-cluster state (see Fig[T]). In 
particular, when the maximum of the sensitivity function is located in the 
second half of the period, the one-cluster (or completely synchronized) state 
acts as a global attractor. In the case, when the sensitivity function reaches 
its maximum in the first half of the period, various two-cluster states become 
stable. 

1.1 Pulse-coupled oscillators 

In some coupled systems, e.g. neuron populations, the time, during which 
the interaction effectively takes place is much smaller than the characteristic 
period of oscillations. In such cases, it is reasonable to approximate the 
interaction by an impact, i.e. by assuming that the interaction is immediate. 
This approximation leads to models of pulse-coupled oscillators, which have 
been widely used in the literature. For example, Mirollo and Strogatz [16] 
have shown, that the complete synchronization (in this case it is equivalent 
to the phase- locking) is stable and attracts almost all initial conditions in the 
network of globally coupled Integrate-and-Fire (IF) oscillators of the form 

dx ■ 

-^ = So-lx,, x,G[0,l), j = l,...,iV (1) 

with constants 5*0 > 7 > 0. One might refer to Sq as input current and to 
7 as the dissipation constant. The following additional condition describes 
the interaction: when k-th oscillator reaches the threshold Xk{t~) = 1, then 
positions of all remaining oscillators are shifted accordingly to the rule 

Xj(t+) = min{xj(t) + X, 1}, j^k (2) 

with some small x > and the A;-th oscillator resets to Xk{t^) = 0. It is shown 
in [16] , that complete synchronization is achieved after a finite transient time. 
The synchronization in a more general model of IF neurons has been shown 
in [T7]. Tsodyks et al. have demonstrated in [18] that the phase- locked state 
is unstable with respect to inhomogeneity in the local frequencies, i.e. when 
the oscillators become nonidentical. 
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Figure 1: Clusters in a population of 50 phase-oscillators. Dots indicate the 
times when an oscillator reaches the threshold; (a) shows the firing pattern of 
a complete in-phase synchronized population (one-cluster), while (b) shows 
the firings in a symmetric two-cluster. 

A larger class of pulse-coupled models was studied in [191 1201 121]. In 
particular, Goel and Ermentrout [12] obtained sufficient conditions for the 
stability of a completely synchronous solution. We introduce this class of 
models in the subsection 11.21 

The dynamics of pulse-coupled oscillators has been studied also for sys- 
tems with different topologies, i.e. ring topology [22], as well as for delayed 
interactions [22] • Transient phenomena of randomly diluted networks have 
been analyzed in [23]. Globally pulse-coupled IF oscillators with a finite 
pulse- width have been considered in, e.g. [251 ISSl |27], where the interaction 
pulse is assumed to have a shape s_*e-"* with the width a. 

1.2 Phase-response curve as a parameter 

In this subsection we introduce a general class of pulse-coupled phase oscilla- 
tors [IHl [28] . The oscillator's motion between the spikes is described by the 
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rule 



(3) 



dt 



where ipj G [0,27r]. When k-th oscillator reaches the threshold at time t, i. 
e. ipk{t~) = 27r, it emits a spike to all other oscillators of the network, which 
are immediately resetted according to 



where Z{ip) is called phase response curve (PRC). Effectively, this means that 
there is no coupling between two consecutive spiking events. The coupling 
occurs only during the spike and acts through the resetting, since the time of 
the resetting of the oscillator j depends on the phase position of the oscillator 
k. The size of the phase-jump, that an oscillator performs, when stimulated 
by an incoming spike depends on its sensitivity to stimulation in its present 
state. See figure [2] for an illustration. 

Let us firstly show, that IF oscillators ([T]) can be written in a form similar 
to ©-(ll]), see also [12]. For this, we rewrite ([T]) with respect to the phase 
coordinate instead of the voltage coordinate. Indeed, the coordinate Xj in 
system ([1]) is supposed to describe the voltage difference across the membrane 
of a neuron [29]. The phase coordinate ipj should behave accordingly to ([3]) 
with the frequency u = 2n/T, where T is the period of oscillations without 
interaction and can be found from ([T]) 




= ^j(t ) + xZ((^j(t )), J ^ k, 



(4) 




The corresponding transformation of variables x = f{^p) can be found from 
the condition 





This gives the function 
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Figure 2: Periodic spiking in a Hodgkin-Huxley neuron model. Solid black 
lines show the evolution of the voltage-component of the model when per- 
turbed by a weak pulse at t = 9 in (a), resp. t = 12 in (b). The dashed 
black lines show the unperturbed oscillations. The PRC Z{if(t)) of the un- 
perturbed model is plotted solid grey and the dashed grey line corresponds 
to Z = 0. The outcome of the perturbing pulse depends on the time t, or 
equivalently on the phase (p{t), of its application. Either, the phase is delayed 
as in (a), i. e. Z{(p(t)) < 0, or it is forwarded as in (b), i. e. Z{ip(t)) > 0. 

which maps the interval < (p < 27t into < x < 1. In the transformed 
coordinates (pj, the dynamics between the spikes is described by (jS]). It 
remains to specify the dynamics at the threshold. Taking into account ([2]), 
when k-th oscillator reaches the threshold ipk{t~) = 2ti its phase resets to 
Vk{t~^) = and all other oscillators have the impact 

^jin = r' {x,{t+)) = /-^(x,(t) + xA(x,x,(t))) 

= rM/((/^,(t)) + xA(x,/(^,(t)))), 

where A(x, x) = min{l, (1 — x)/x} < 1. In the case of small x, i.e. the 
assumption of weak coupling holds, the resetting rule can be approximated 
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Figure 3: Phase- response curve for IF model ([1]). The function Zip^Lp) 
measures the sensitivity of the system to a small external perturbation at 
different positions ip. The corrected function Zjp^^^ip) does not allow the 
oscillators to be moved over the threshold by a spike. 

as 

V9j(t+) = ^pj{t) + xmm.{ZiF{^j{t)), (27r - ifj{t))/x}, (6) 

where 

Z^^) := -^im) = ^exp ^-^j . (7) 

Thus, with respect to the phase coordinates, the IF model ([T]) has the form 
([3]), (E]). In particular, the resetting rule is given by the function 

ZifA'^) = niin{Z/i7'((^j(t)), (27r - (fj{t))/x}, (8) 

which depends on the amplitude of the perturbation x. Figure [3] illustrates 
this function for x = 0.05. Practically, the PRC measures the sensitivity of 
the phase to external perturbations. 

We have shown above the specific example of pulse-coupled IF models 
and their reduction to pulse-coupled phase oscillators Q-dl]). In fact, this 
procedure is also possible for higher- dimensional smooth systems, whenever 
the oscillations correspond to a hyperbolic limit cycle, i.e. in a generic case. 
More details can be found in [191 1211 EO]. When the coupling is acting along 
one component, e.g. the voltage variable, as often assumed in the case of 
neural populations, the PRC appears as a scalar function of the phase. In 
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Figure 4: Examples of different PRCs, (a) Hodgkin-Huxley model, (b) 
Connor model. Note that the functions and their derivatives are zero at the 
ends of the interval 99 = and ip = 2ti (adapted from [SI]). 



the case of a higher-dimensional interaction, it should be considered more 
generally vector. 

Examples of PRCs for different neuron models are shown in Fig. HI Some 
more numerically and experimentally obtained PRCs can be found in e.g. 
jl9| [28l I3T] . The remarkable feature of many of such PRCs is that, contrary 
to the IF model, their PRCs are independent on x and admit zero values 
at (/9 = and ip = 2tx. The conditions ^(0) = Z(27r) = are also reason- 
able from the neuroscientific point of view, since they reflect the fact that the 
neurons are not sensitive to perturbations during the spike (see Fig. H]). Gen- 
erally speaking, system (I3])-(I1]) is a useful model, which possesses quite a big 
generality by including the PRC as some "infinite- dimensional" parameter. 



1.3 System description 

Our main object of study is the following system of globally pulse-coupled 
phase oscillators of the form 

rlin ■ 

(9) 



dipj 



dt 

with the resetting rule 

Mt-^) = 0; <^i(t+) = ^jin + ^^(^,(r)), J ^ k, (10) 
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Figure 5: Family of the unimodal PRCs see f|TT]) . 

where the velocity of the phase is assumed to be 1 without loss of generality. 
We assume a fixed, positive overall coupling strength x > 0. The impact is 
rescaled taking into account the number of oscillators, see also [27]. In this 
study, we consider a one-parametric family of the PRCs, which are positive 
and unimodal as shown in Fig. [51 The parameter /3 G [0, 1] controls the 
position of the maximum, namely, for larger /3, the maximum is located in the 
domain of small </:>, which corresponds to a more sensitive excitatory response 
of the system just after spike. For smaller /3, the system is more sensitive 
to perturbations shortly before the spike. The value /3 = 0.5 corresponds to 
an intermediate situation. We assume also that Z'{0) = Z'(27r) = 0, which 
is appropriate for a broad class of experimental and analytically obtained 
PRCs (see Fig.SD. 

We note that the qualitative results reported in the paper are independent 
on the exact expression for the PRC but rather on the shape of the PRC and 
its behavior aX (f = and (p = 2ti. Our particular choice is 

Z^(v.) = l-cos^^(y.), /3e[0,l], (11) 

where 

In particular, Z^^^^if) = 1 — cosy^. 

System (l9l)-(fT0l) is equivalent to an (A^ — 1) -dimensional discrete dynam- 
ical system, which can be obtained as a return map by considering its state 
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each time when some of the phases reaches a fixed value, e.g. = 27r. Let us 
point out how this map appears. Without loss of generality, we may assume 
that the phases are ordered as 

271 = (fi> ip2> ■ ■ ■ > (12) 

at t = 0. We will use the important property of fl9l)- f|T0|) that the oscillators 
do not overrun each other for all times if the system size N is sufficiently 
large. Indeed, since the inequality 

+ J^^^^i) ^ Vj+i + j^Z (v^j+i) (13) 

holds for sufficiently large A^, the order of oscillators is preserved during the 
spike. It is also evident, that the order is preserved between the spikes as 
well. More exactly, the inequality 2ti > ipi^i > (p2+i > ■ ■ ■ > (pN+i > holds 
for all t, where / is some shift and the indices are considered modulo A^. 

Let us denote by Ki the map, which maps the initial phases ( lT2l) into 
the phases at the moment when the oscillator (p2 reaches the threshold, i.e. 
ip2 = 27r. It is easy to obtain that 



Ki{(fi,(f2,f3, ■ ■ • ,V57v) = 

= (27r - /i(v32), 27r, /i(v93) + 2%- /i((/?2), • • • , /^(v^jv) + 2n - fi{ip2)), 

where 

K^):=V + ^Ziip). (14) 

In a similar way, the mapping K2 exists, which maps the phases to the state, 
where the third oscillator is at the threshold and so on. The composition of 
maps 

K = Kno Kn-i o---oKi (15) 
gives the dynamical system on A^- dimensional torus 

(V^i, . . . ,V3Ar) i^(V9i, . . . ,V?7V), (16) 

which maps the initial state fll2p into a new state after all oscillators have 
crossed the threshold once and the first oscillator reaches again the threshold. 
We call the map K return map. 
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In this paper, we will not use the explicit form of the mapping f llSp . For 
our purposes it is important to conclude that the dynamics of system ([H])- 
ffTOl) are indeed equivalent to some (A^ — l)-dimensional, discrete dynamical 
system on the iV- dimensional torus. The smoothness of this system depends 
on the smoothness of its PRC function. 

2 Numerical results 

In order to detect the appearance of one- or two-cluster states, we have 
numerically computed the order parameters 



A perfect one-cluster state is characterized hj Ri = R2 = 1 and a perfect 
antiphase two-cluster is characterized by i?i = and R2 = I. We present re- 
sults of simulations for x = 0.5, but qualitatively we observe similar behavior 
for a broad range of x > 0. 

As shown in Fig. El we observe two qualitatively different types of be- 
havior depending on parameter /3. For /3 < 0.5, i.e. when the maximum of 
the PRC is shifted to the right (see Fig. [5]), the one-cluster state seems to 
be the attractor; for /3 > 0.5 and the maximum of the PRC is shifted to the 
left, a two-cluster state is attracting. We have chosen initial conditions in a 
vicinity of a two-cluster state in Fig. El^a) and (b), therefore the initial values 
of the order parameters are i?i ~ and i?2 ~ 1. Figure [6](b) shows how the 
instability of the two-cluster state implies desynchronization transient, after 
which the system is attracted to a synchronous one-cluster state. Similar 
behavior occurs for other initial conditions. Figure El^c) and (d) illustrate 
the order parameters behavior for initial conditions close to the splay state 
(a state, where the phases are distributed). The initial values for the order 
parameters in the splay state are close to zero, but after a transient, they 
approach again the same asymptotic values as in (a) and (b). 

A more complicated behavior occurs for the intermediate value of the 
parameter /3 = 0.5, i.e. when the PRC is symmetric. In this case, the order 
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Figure 6: Behavior of the order parameters Ri(t) and R2(t) for a trajectory 
starting in a vicinity of the two-cluster state for (a) and (b). The lower panel 
(c) and (d) corresponds to a trajectory starting in a vicinity of the splay 
state. Left figures (a) and (c) correspond to the parameter value /3 = 0.7, 
where the two-cluster state is attracting and (b) and (d) to /3 = 0.3, where 
one-cluster state is attracting. 



11 




(b)l 



0.75 



1 1 


1 _ 

1 1 





0.5 


1 



Figure 7: Dependence of the asymptotic values for the order parameter Ri 
(a) and R2 (b) on (3. For the most values of /3, except (3 = 0.5, the order 
parameters tend to some constant value, when initialized near the splay state 
or the symmetric two-cluster state. 



parameters Ri{t) and R2{t) do not approach some asymptotic constant values 
but remain periodic in time. As a result, the maximum asymptotic values 
of both Ri and R2 do not coincide with the corresponding minimum values. 
This type of behavior is observed for a very small parameter interval of order 
10~^ around f3 = 0.5. We discuss it in Sec. |5] in more details. Figure [7] 
summarizes the behavior of the order parameters for different /3. 



3 Appearance and stability properties of one- 
cluster state 

In an ideal one-cluster synchronized state, all oscillators have the same phases 
(fj = (fs for all j. This state is a fixed point of the map (1161) . because the 
PRC turns to zero at = 27r and (f = 0. This means that the coupling 
vanishes for one-cluster state. More exactly, when an oscillator ipj fires, i.e. 
ipj = 27r, all other oscillators have the phase 27r and do not obtain the spike. 
As a result, the period of this state is determined simply by the uncoupled 
dynamics and equals 2%. 
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3.1 Inadequacy of the linear stability analysis 

In order to obtain conditions for the stability of one-cluster state, one can 
examine the return map f|T6|) . The linearization of this return map around the 
one-cluster state gives then the corresponding multipliers, which determine 
its local linear stability. As it is expected, the local stability is governed by 
the properties of the PRC a.t if = and if = 27t. This procedure has been 
done in [19]. Applying these results to our case, the resulting conditions for 
the local linear stability of one-cluster state is 

1 + -Z'i2n-)) (l + -Z'iO^)) <1, / = 1,A^-1. (19) 

We observe that the necessary condition for the linear stability is that the 
derivatives of Z{ip) at the ends of the interval [0, 27r] do not vanish. This is 
not the case for our PRC fITT]) . Hence, all associated multipliers have modulus 
one and the linear stability analysis do not provide useful information about 
the stability of one-cluster state. 

3.2 One-cluster state is a saddle point 

In this section we show that one-cluster state is a saddle point, i.e. there are 
some arbitrary small perturbations of this state, which grow with time. At 
the same time, some other small perturbations decay. 

Existence of a local unstable direction. First of all, let us show that 
one-cluster state is unstable with respect to the following special perturba- 
tion: 

(fi = (fs + e, ^2 = ■ ■ ■ = = (20) 

with arbitrary small e > 0. During the period between spikes, the dynamics 
is monotonous ^j{t) = ip^ + t ioi j = 2, . . . , N and (pi(t) = ips + £ + t, thus, 
the distance between the phases remain constant. Without loss of generality 
we may assume that 

V9i(0-) = 2n, ^2{0-) = ■■■ = (^Ar(O-) =27i-e. 

After the first oscillator moves over the threshold and resetting occurs, the 
phases follows 

(^i(0+) = 0, (^2(0+) = ■ ■ ■ = <^^(0+) = 27r-e + ^Z{2Ti - e) = ^(27r - e). 
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The next resetting occurs at time ti = 27r — 7^2(0^) = e — j^Z(27r — e) when 
the group of — 1 synchronous oscillators reaches the threshold. At this 
moment 

¥>i{tl) = e - ^Z{2tx - e) > 0, (^2(^7) = " " " = ^a^(V) = 2vr. 

Now the group of — 1 synchronous oscillators is at the threshold. The 
correct definition of the firing rule for this case can be naturally obtained by 
extending it to the situation when all the oscillators V92, • . . , v^at in the cluster 
have slightly different phases and then allowing the phases to converge to 
the same value. This leads to the following resetting rule when passing the 
threshold by the A^ — 1 cluster: 

^M) = /x^-^ (^i(tr)) = /i^-^ {e - ^Z(2vr - e)) , (21) 

^^{tt) = --- = ^N{tt) = Q, (22) 

where denotes the superposition of A^ — 1 functions /io/io/io---o/i, 

where is defined by f|T^ . The resetting fl2T|) simply means that the function 
/i is applied A^ — 1 times (whenever an oscillator from the cluster Lp2i ■ ■ ■ ,fN 
fires) in order to obtain the final position of (fi. 

In this way, we obtain a mapping, which maps the initial size of the 
perturbation e at time t = into its new size Yi{e) at time ti. The mapping 
is 

e ^ Y,ie) = /i^-i [e - ^Z{2n - e)) . (23) 

It is clear that ^1(0) = 0, what corresponds to the invariance of the one- 
cluster, and the stability properties of the origin of ( 1231) determine the sta- 
bility of the one-cluster state with respect to the specific perturbation fl20!) 
chosen. Up to the linear level, the origin of fl23|) is neutrally stable, i.e. 
Y({0) = 1, which is clear, since the one-cluster state is linearly neutrally 
stable. The second derivative of (1251) at e = is nontrivial 

Y;'{0) = xZ"(0) - ^ (Z"(2vr) + Z"m 

and is positive for sufficiently large A^ since Z"{0) > for /3 G (0, 1]. Hence, 
for sufficiently large A^, the origin of fl23l) is unstable, see Fig.|8l^a). This leads 
to the local instability of one-cluster state for all (5 G (0, 1]. Accordingly to 
this, the distance of the advanced oscillator (pi from the remaining cluster 
will grow, but this growth is not exponential. 
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Existence of a local stable direction. Now let us show that the one- 
cluster state is locally stable with respect to perturbations of the form 

Lpi = LPs- e, ^2 = ■ ■ ■ = = (24) 

with £ > 0. This can be shown similarly to the previous case by obtaining 
the discrete mapping, which describes the dynamics of the perturbation. In 
the case of perturbations flMl) . this mapping reads 

e^y7v-i(£) =yu(27r-/i^^^(27r-£)) (25) 

and has the following properties 

>Ar-l(0) =0, 

and 

r;;_,(0) = -xZ"(27r) + ^ iZ"{2n) + Z"(0)) (26) 

It implies that for sufficiently large N the second derivative is negative and 
the origin of the discrete mapping e — )■ Ypi^i^e) is asymptotically stable (see 
Fig. [8](b)). Hence, the one-cluster state is stable with respect to perturba- 
tions of the form (!24|) . This, together with the instability with respect to 
perturbations ( |20l) . implies that the one-cluster state is the saddle point in 
the phase space (see schematically Fig. [9]). 

Other stable and unstable local directions. In general, the two-cluster 
perturbations of the one-cluster state are given by 

if,{0) = ... = <^^j0)=27r (27) 
<^^,+i(0) = ... = <^^(0) = 27r-£, (28) 

where A^^i -|- A'"2 = A^. This means, there are A^^i oscillators in the front- 
group and the remaining N2 oscillators in the back-group. The corresponding 
discrete 1-D systems, which describe the dynamics of such perturbations are 
given by 

e — )■ YN^{e) and e YN^{e), 
where Yj{0) = 0, ^Fj(O) = 1 for j = 1, - 1 and 

^^YnM = ^ {N2Z'\0) - N,Z'\2n)) , (29) 
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Figure 8: Local Cobweb-Diagram of the functions Yi{e) and Yi^_i{e) around 
e = 0. Iterations of these maps determine the behavior of special perturba- 
tions to the one-cluster state, (a): small perturbations grow with time; (b): 
small perturbations decay. 

^2YnM = ^ (iViZ"(0) - N,Z"{2n)) . (30) 

The expressions ( 129|) and ( 130|) may have different signs depending on the 
values of Ni, N2, as well as the second derivatives Z"{0) and Z"{2tt). This 
implies the existence of multiple unstable as well as stable directions to the 
one-cluster solution, for more details, see section HI 

3.3 Stable homoclinic orbit to one-cluster state 

Let us first note that the two-clusters of the form f E7|l - fl25]) do not split with 
time. In geometric terms, this means, that the subspace corresponding to 
such solutions is invariant. In particular, the subspace, which corresponds to 
A'"! = 1 and A''2 = iV — 1 is invariant as well. Being restricted to this invariant 
subspace, the one-cluster state is a saddle point, as we have shown in the 
previous section. In Appendix [7] we prove that there exists a homoclinic orbit 
in this subspace, which connects the both unstable and stable manifolds, see 
Fig. [9l In fact, as will be shown in Sec. HI the dynamics within the invariant 
subspace is given by the 1-D mapping shown in Fig. [TlT b). 

Numerical calculations further supports this result and show that the 
invariant set, which is composed of a homoclinic loop and the fixed point is 
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Figure 9: One-cluster state as a saddle point in the phase space with a 
homoclinic loop. 

an attractor. Figure HD] shows how the width of the cluster changes as time 
evolves for some typical initial conditions. More specifically, we compute 



One can clearly observe that the width tends eventually to zero interrupted 
by some blowouts. The blowouts correspond to the events, during which the 
first oscillator leaves behind the remaining cluster and makes a rotation in 
the phase. After the rotation, it joins again the cluster and becomes the 
"last" one. The time interval between such events grows unboundedly with 
time supporting the homoclinic nature of the attractor. Note that the width 
of the cluster should be nonzero in order to observe this phenomenon, i.e. 
one should perturb the system slightly from the fixed point, see Fig. [91 

Finally, we would like to remark that the same methods allow proving 
the existence of other homoclinic orbits, which correspond to two-cluster 
perturbations (ITTjl - fl^Sl) with Ni <^ Hence, one should rather speak 
about an attracting family of homoclinic orbits. 

4 Two-cluster states 

Two-cluster state appears when the oscillators split into two groups (see 



A{t) 



max 



Fig. [ID 



(31) 
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Figure 10: Width of the cluster A(t) = maxi<ij<N {\ipi(t) — ipj(t)\} as a 
function of time. Figure (a) shows the behavior along the orbit started at an 
initial condition close to the splay state (far from the one-cluster), (b) shows 
the behavior along the orbit started close to the state f l20l) . The behavior 
indicates the existence of a stable homoclinic orbit. 
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Contrary to one-cluster state, the two-cluster state must not be a fixed 
point of the return map f ll6p . Indeed, when two clusters appear, their relative 
behavior is then given by the following discrete return map (by assuming that 
the return map is computed for ipi = 2tt and ip2 < i^i) 

ip2 ^ Yn,{^2) ■■= 27r - /i^^ (27r - /i^H^2)) • (32) 

This map has different properties depending on Ni, N2 = N — Ni as well as 
on /3. All such maps have zero fixed point corresponding to the case when 
two clusters merge into one. One can obtain 

rjvi(0) = 0, F;vi(27r) = 27r, 

r^(0) = l, F^(27r) = l, 

and 

Y^2n) = ^{N2Z"{27r)-N,Z"m. 

Figure [TT] shows typical maps for three different situations: 

(a) The map has an unstable fixed point inside the interval [0, 27r] and the 
endpoints x = and x = 27r are asymptotically stable. Hence, within the 
corresponding subspace, the one-cluster state is asymptotically stable (simi- 
larly to Fig. El^b)). 

(b) The map has unstable fixed point at x = and stable at x = 27r. This 
case corresponds exactly to the case, when the one-cluster state has a homo- 
clinic orbit starting in x = and ending at x = 27r (0 ~ 27t on the torus). 

(c) The map has a stable fixed point inside the interval [0, 2tt] and the end- 
points X = and x = 27r are unstable. Hence, within the corresponding 
subspace, the one-cluster state is asymptotically unstable and the two-cluster 
stationary state is stable. 

The fixed points of the map fl32|) give two-cluster stationary states: 

i^2 = YN,M- (33) 

The condition for the merging of two cluster into one cluster is given by 
the condition for the existence of the double root of the function Yn-^{iIj) at 
^ = or ^ = 27r, i.e. Yl^^{0) = or F^^(27r) = 0. This results into 

NiZ"{0) = N2Z"{2n). (34) 
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X 271 X 2n x 2n 



Figure 11: Typical behavior of functions Yfyr^{x) (see fl32l) ). which determine 
the behavior of two-clusters. On the figure Ni = 150 and = 500. 

Expression ( l34l) determines also the moments when one-cluster state un- 
dergoes pitchfork bifurcations. At such bifurcation, two different nonsym- 
metric two-clusters bifurcate from the one-cluster state: one with A'^i = pN, 
N2 = {1 — p)N, and another with A''i = {1 — p)N, N2 = pN . The bifurca- 
tion diagram in Fig. [12] shows some of the branches of two-clusters, which 
originate from %l)2 = Q 01 11)2 = 27r. 

The pitchfork bifurcations for /3 < 0.5 are subcritical. Namely, the two- 
cluster states are unstable and they merge into the one-cluster state. With 
increasing /3 the one-cluster state becomes more and more locally unstable 
by transforming stable directions into homoclinics (see Fig JTTl) . In spite of 
this fact, we observe numerically, that the invariant set, which is composed 
of the one-cluster state and homoclinic connections is still attracting in the 
phase space. All two-cluster states, which exist at this moment, are unstable. 
As a result, one computes high values of the order parameters i?i and R2 on 
the numerically obtained figure [7] for (3 < 0.5. 

4.1 Stability of two-cluster states 

For P > 0.5, the invariant set composed of one-cluster state and homoclinic 
orbits losses its stability and two-cluster states emerge, which are asymptot- 
ically stable. Numerical results in Fig. [13] show which two-clusters are stable 
depending on the parameter /3. In general, for P closer to 0.5, the symmet- 
ric clusters with p ^ 0.5 are stable. As /3 increases, the more asymmetric 
clusters stabilize as well. This implies that the PRCs with the maximum. 
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■p=0.51 
p=0.7 

■p=0.98 




P 



(b) 

271 



-p=0.498 
-p=0.5 
p=0.502 



/3 



Figure 12: Positions of the two-cluster states 6 = 2n — ip2^ where il)2 are fixed 
points of f l33l) . Different fines correspond to different cluster splittings, i.e. 
A'^i = pN , N2 = (1 —p)N. At (5 = or (5 = 27r, tfie corresponding two-cluster 
is merging into tfie one-cluster. 



wfiicfi is sfiifted to tfie left favor tfie coexistence of a large number of stable 
brancfies of two-clusters. 



5 Intermediate state for symmetric PRC with 

p = 0.5 

Tfie case of symmetric PRC for /3 = 0.5 is degenerate. Wfien increasing 
/3 tfirougfi 0.5, tfie fiomoclinic sets including tfie one-cluster state become 
unstable and a two-cluster state becomes stable as it is described in tfie 
previous section. Tfie numerical calculations for /3 = 0.5 sfiow nonstationary 
dependence of tfie order parameters on time, see Fig. [TH One observes 
periods of time, wfien two-clusters persist. Tfiese periods are cfiaracterized 
by almost constant order parameters. Tfie periodic blowouts of tfie order 
parameters correspond to tfie befiavior, during wfiicfi tfie oscillators from 
tfie advancing cluster spread over a big part of tfie pfiase circle and finally 
form anotfier cluster befiind (see tfie inset in Fig. [T^ . 
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(b) 




Figure 13: Stability and existence of two-cluster states, (a) Solid lines de- 
note stable two-cluster stationary states and dashed - unstable. The lines are 
shown only for selected values of p = Ni/N, while the dense set of branches 
for all possible p exist. Figure (b) shows which two-clusters are stable in de- 
pendence on /3 (obtained numerically), p = 0.5 corresponds to the symmetric 
cluster and p ^ 0.5 to nonsymmetric clusters. 




5.000 



Figure 14: Nonstationary behavior of the order parameters -Ri and R2 with 
time for /3 = 0.5001. One observes periodic restructuring of two-clusters. 
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6 Conclusions 



In this paper we have studied the asymptotic behavior of a system of globally 
pulse-coupled phase oscillators (!9|)- (fT0|) with the phase response function, 
which is positive, unimodal, and turns zero at the threshold together with its 
first derivative. In particular, we considered the question how the position 
of the maximum of the PRC influences the dynamics of the coupled system. 

We have numerically observed, that for the PRCs with the maximum 
shifted to the right (for our model, it corresponds to /3 < 0.5), a one-cluster 
state becomes apparently stable. More detailed analysis reveals that the 
one-cluster state is, in fact, asymptotically locally unstable, i.e. a generic 
small perturbation will grow with time. Moreover, we show that trajectories 
of the system has a behavior, which is characterized by long-time intervals 
when the system stays close to the one-cluster state and long excursions 
away from the one-cluster state (see Fig. [TOj) . The excursions become less 
and less frequent with time. This behavior is explained by the existence 
of the family of homoclinic orbits to the one-cluster state, which altogether 
form an attracting set in the phase space of the system. 

In the case, when the maximum of the PRC is shifted to the left, i.e. 
the oscillators are mostly sensitive to perturbations in the phase just after 
the threshold, the one-cluster state no more dominates the dynamics and 
various stationary two-cluster states become stable. These two-cluster states 
appear in pitchfork bifurcations from the one-cluster state as parameter /3 
increases. First, at /3 = 0.5, there appears a symmetric two-cluster with equal 
number of oscillators in each cluster. With further increasing /3 more and 
more asymmetric clusters appear and become stable leading to the increasing 
coexistence of stable two-clusters. 

7 Appendix: Existence of a homoclinic orbit 



Theorem. For (5 G (0, 1) there exists Nq, such that for populations of size 
N > Nq, system / f75]) possesses a homoclinic trajectory, which connects the 
one-cluster stationary state. The homoclinic trajectory has the form 



where lim. 



-n— >■— oo 



(pi{n) = 0^ and lim„^+oo </5i('^) = 27r~. 



(35) 



23 



Proof. Fix P G (0, 1) . We will consider 

Yi {N, x) = 27r - /i (A^, 2tt - /i^^^ (A^, x)) , 
where fi^ {N, x) denotes the j-th iteration of 

X h-^ n{N,x) = X + —Zp (x) . 

The map Yi {N,x) describes the evolution of the distance x G (0,27r) during 
a time interval in which all oscillators of a population (pi = x; ip2 = ■■■ = 
(fN = 27r, emit exactly one spike. Homoclinicity then is equivalent to 

Y^ {N, x) 2tt, for all x G (0, 27r) , as oo, 

where Y^ {N, x) denotes the k — th iteration of x i— )■ Yi (A^, x) . Analogously 
to the analysis of section HJ we find that 

ri(A^,0) = 0; Fi(Ar,27r) = 27r, 

Y;{N,0)=Y;iN,27r) = l, 

YI' {N, 0) > 0, Yl' (N, 2it) > 0. 

Here and in the following, primes denote the derivatives with respect to the 
second argument (phase). For fixed A^, there exists a rejecting region (0,£:Ar) 
where Y(' {N,x) > 0, for x G (0,6^) and an attracting region (27r — 6^,271) 
with YI' {N, x) > 0, for a; G (27r - en, 2n) . This gives: 

Y.'^ {N,x)>eN, 

for X G (0, En) and some finite ko = ko {N, x) G N, and 

Yl" {N, x) 2tt, 

for /c — )■ cxD and x G (27r — En, 27r) . Our goal is to show, that there exists a 
uniform > 0, such that for all N > Nq : 

Yl' {N, x)>0 for X G (0, ^o) and 
Yl' {N, x)>0 for X G {eq, 27r - Eq) , 

and such that for all N > Nq and all x G [eo,271 — eq] : 

Fi(A^,x) > x + Ajv, 



24 



with 



'■= min Yi {N, x) — x > 0. 

a;e[ejv,27r— £]v] 



Thus, any x G (0, 27r) will reach the attracting region (27r — Eq, 27r) within a 
finite number of iterations of x H- Yi (A^, x) . Let us write 

Y,{N,x) = Yi{x) + ^w{N,x), 

where Yi (x) = x + [x) is independent of A^. For Yi we have 



This implies for 



that 



Yx [x] > X for X e (0, 2it) , 
Fi (0) = 0, Fi (27r) = 27r, 
F;(0)=y'/(27r) = l, 

w {N, x) = N (Yi (iV, x) - Yi (x)) 



w {N, 0)=w {N, Itx) = 0, 
w' {N, 0) = w' {N, 27r) = 0. 

We will show, that the region [0,£:7v] may be chosen as [0,£o] i independently 
on large N. The analysis for the other region [27i — Eq, 27r] can be done simi- 
larly. Around x = 0, we have the following representation of Yi (A^, x): 

2 

Y, {N, x) = n {N, 0) + YI {N, 0)x+ ^Yl' {N, 
= Y, (0) + Y; (0) X + ^Yl' (e^) + j^(w {N, 0) + w' {N, 0) x + ^w" {N, ^iv)) 

for some G [0,£:] . Further it holds Y{' (0) > 0. This means, there exists 
an Eq > 0, such that for x G [0,5o] ? ^i" i^) > 0. Now we construct an A^- 
independent lower bound for w" {N, x) in x G [0, Eq] , where Eq will be further 
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altered in the analysis without always choosing a new notation. In other 
words, we claim that there exists Cq G M with 



liminf ( min w" (N,x) ) > Cn. 

N^oo VxG[0,eo] 



(36) 



We have 



w {N, e) = N (Yi {N, x) - Yi (x) 
= N {27r- fi {2n - {N, x)) - x - xZ^ (x)) 
= N (/i^-i (iV, x) - ^Zp (27r - ^^"1 (iV, x))-x- xZ^ (x)) 



Ar-2 



= X V (Z^ [li^ {N, x)) - Zfi (x)) - KZf, (27r - /i^'^ (iV, x)) - xZ^ (x) . 

.=0 ^ 2 ' 

Since part I, as well as its derivatives, is obviously uniformly bounded in 
and X, we restrict us to establish (136|) for 



Af-2 



W 



(iV, x) = Y, [Zp ifi' {N, x)) - Zp (x)] . 
i=o 



We have 



N-2 



w' (AT, x) = ^ [z^ (/i^- (iv, x)) (/i^- (iv, x))' - z;, 



,x 



j=0 



N~2 



w"{N,x) = 5^[z;;(/x^(iv,x))((/iMiv,x))' 



(37) 



j=0 



+Z',{^^^N,x)) {^^HN,x))"-Z;ix) 



To handle this, we need some uniformity-properties of fi^ {N, x) . Elementary 
calculations give 



(/.^- (iv, x)) ' = n {N, {N, x)) = n (i + ^4 (/ ^)) 



fc=0 



fc=0 
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i i-i 



1=0 k=0,k=il 

x^Z'^{^^^{N,x)) (/i'(iV,x))'. 
This implies that the following inequality 

< {N, x))' < exp (xCO , where C' = max \Z'. (x) I (38) 

a;G[0,27r] 

holds for all large enough A^. This again yields 

X < fi^ (N, x) = /i^iV, 0) + /" {fi^ (AT, y)ydy 

Jo 

< x + xexp{KC'). (39) 

Using this upper bound, we get some A^- independent Eq, such that for x e 
[0,eo] 

Z'^ifi'' (iV,x)) >0. 
This gives A^-independent monotonicity of 



X 



^ Z'^ [fi'' {N, x)) > for X e (0, Eo) . 



Further, we can use 0391) to improve the bounds fl38p for (/i-' (N^x))' in x G 
[0,£o] to 

1 < (/i-'' (A^, x))' < exp (xC') . (40) 

This implies 

/i-' (A^, x) < X ■ exp (x^') . 

We find 

< {fi^ {N, x))" < exp (2xC') < xC" exp (2xC') , 



where 



Now observe that 



C" = max (x) . 

xe[0,27r] ' ^ ^ ' 



Z';{0) + Z'^ (0) = {A - -) + -(3 + - > 0, 

Tf J IT IT 
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i.e., eventually further decreasing of eo > gives, with Eq = Eq ■ exp {xC) '■ 

(41) 



min Z'f^ (y) > - mm Z'^ (y) . 

yG[0,£o] ye[o,eo] 



Hence, for a; G [0,Eo] : 



Af-2 
j=0 



>0 



i=o 

N-2 



i=o 

where we have omitted the arguments (A^, x) of fi for brevity. Using (HTI) . we 
continue the estimations 

/ / 2 \ \ 



Af-2 






2 


( 




\ 


min Z« (y) > 


(1 




-1- r 

Jo 




-1 


dy 








) 


I 



N~2 



j=0 



N-2 



> min Z;;(y)^(l-a:)((/i^y-l) >0. 



j=0 



This establishes ( 136|) and hence Yi (A^, x) > x G [0, £o] for large enough N. □ 
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